function [ skewerStats ] = getSkewerStats( sktp, length, rowsInFrame )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
    
 
      
   %% Transform subscripts to positive x and y axis (upper left quadrant).
   
   maxRow = rowsInFrame;
   coords = zeros(size(sktp));
   coords(:,1) = sktp(:,2);
   coords(:,2) = maxRow - sktp(:,1);
   coords = transpose(coords);
   numPoints = size(coords,2);

   %% Calculate the rotation angle 
    xHead = coords(1,1);
    yHead = coords(2,1);
    xTail = coords(1, numPoints);
    yTail = coords(2, numPoints);
    yDiff = yHead - yTail;
    xDiff = xHead - xTail;
    
    % ...head is pointed in upper left 
    if xDiff < 0 && yDiff >= 0
        
        theta = - atan( abs(yDiff/xDiff));
        skewerAngle = 180 - atand(abs(yDiff/xDiff));
    % ... head is pointed to the upper right
    elseif xDiff > 0 && yDiff >= 0 
        theta = -(pi - atan( abs(yDiff/xDiff) ));
        skewerAngle =  atand( abs(yDiff/xDiff) );
     % ... head is pointed to the botton left
    elseif xDiff < 0 && yDiff <= 0 
        theta = atan( abs( yDiff/xDiff)) ;
        skewerAngle = 180  + atand( abs(yDiff/xDiff) );
     % ... head is pointed to the botton right
    elseif xDiff > 0 && yDiff <= 0 
         theta = - atan( abs(yDiff/xDiff) );
         skewerAngle = 360 - atand( abs(yDiff/xDiff) );
    end
   
    R  = [cos(theta) -sin(theta)  ;...
        sin(theta)  cos(theta)   ];
        
    %%  Analysis with centroid at the origin
    
   
    
    rotatedCoords = R * coords;
    
    %Shift the skeleton so the centroid is at (0,0)
    xCentroid = sum(rotatedCoords(1,:) )/numPoints;
    yCentroid = sum(rotatedCoords(2,:) ) /numPoints;
    transformedCoords = rotatedCoords -repmat([xCentroid;yCentroid], 1, numPoints);
    
   
   
    
   %Get the min.max x/y coords
    xMax = max(transformedCoords(1,:));
    xMin = min(transformedCoords(1,:));
    yMin = min(transformedCoords(2,:));
    yMax =  max(transformedCoords(2,:));  
    
    %Get distances from the x axis to the largest and smallest y values.
    A = abs(yMax);
    B = abs(yMin);

    sktHeight = A + B;
    sktWidth = abs( xMax - xMin);
    sktAmpRatio = min(A,B)/max(A,B);
    sktCmptFactor = length/( sktHeight* sktWidth);
    sktElgFactor = sktHeight/sktWidth;

      
    Ixx = sum(transformedCoords(1,:).^2)./numPoints;
    Iyy = sum(transformedCoords(2,:).^2)./numPoints;
    
    %Get the average radial angle from each skeleton point to the centroid
    %using element-wise arithmetic
    phi = atand(transformedCoords(1,:) ./transformedCoords(2,:));
    aglAve = sum(phi)/numPoints;
    
    
    %% Analysis with origin the point 1/2 dist. between the min and max row and col
    
    %Rotate the points
    rotatedCoords = R * coords;
   
    xCenter = (  min(rotatedCoords(1,:)) + max(rotatedCoords(1,:))  )/2;
    yCenter = (  min(rotatedCoords(2,:)) + max(rotatedCoords(2,:))  )/2;
    
    %nsform the skeleton so the midpoint is at (0,0)
    transformedCoords2 = rotatedCoords -repmat([xCenter;yCenter], 1, numPoints);
    
    
    xsym = sum(transformedCoords2(1,:));
    ysym = sum(transformedCoords2(2,:)); 
    xysym = sum( (transformedCoords2(1,:).* abs(transformedCoords2(2,:))));
    
    trackAmplitude =  max( transformedCoords2(2,:)) - min(transformedCoords2(2,:)) ;
    
     
    %Frequency Analysis to get characteristic period
     
    signal = transformedCoords(2,:);
    [maxValue,indexMax] = max(abs(fft(signal-mean(signal))));
    
    %Length of signal
    L = numPoints;
    frequency = indexMax / L;
    period = 1/frequency;
    
    
    
    
    
    %% Load the results
    
    skewerStats.SktHeight = sktHeight;
    skewerStats.SktWidth = sktWidth;
    skewerStats.SktAmpRatio = sktAmpRatio;
    skewerStats.SktCmptFactor = sktCmptFactor;
    skewerStats.SktElgFactor = sktElgFactor;
    skewerStats.Ixx = Ixx;
    skewerStats.Iyy = Iyy;
    skewerStats.SktAglAve = aglAve;
    skewerStats.Xsym = xsym;
    skewerStats.Ysym = ysym;
    skewerStats.XYsym = xysym;
    skewerStats.TrackAmplitude = trackAmplitude;
    skewerStats.TrackPeriod = period;
    skewerStats.SkewerAngle = skewerAngle ;
    
       
end

